function KF0(H_t,Se,Phi1,y,n_agg)

(T,n) = size(y)
s     = zeros(T,n);
P     = zeros(T,n,n);

# Initialization
s[1,:]   = y[1,:];
P[1,:,:] = H_t[:,:,1];;

# Kalman Filter Recursion
sprime   = zeros(n,1);
liki     = zeros(T,1);

# Forward Iterations
for t=2:T

    # Updating Step
    sprime = Phi1*s[t-1,:];
    Pprime = Phi1*P[t-1,:,:]*Phi1' + Se;

    # prediction step
    yprediction = sprime;
    v = y[t,:] - yprediction;
    F = Pprime + H_t[:,:,t];
    # if det(F) <= 0.0
    #     U,Sig,V = svd(F)
    #     detF = prod(Sig[Sig.>1e-12])
    # else
    #     detF = det(F)
    # end
    liki[t]    = -0.5*n*log(2*pi) - 0.5*log(det(F)) - 0.5*v'*inv(F)*v;

    # updating step
    kgain    = Pprime*inv(F);
    s[t,:]   = sprime + kgain*v;
    P[t,:,:] = Pprime - kgain*Pprime;

end

return s, liki

end

#-----------------------------------------------------------

#-----------------------------------------------------------

function KFnoX(H_t,Se,Phi1,y,n_agg)

(T,n) = size(y)
n_s   = n-n_agg
s     = zeros(T,n_s);
P     = zeros(T,n_s,n_s);

y_obs = y[:,1:n_agg]
s_obs = y[:,n_agg+1:end]

# partition PHI1
Phi_yy = Phi1[1:n_agg,1:n_agg]
Phi_ys = Phi1[1:n_agg,n_agg+1:end]
Phi_sy = Phi1[n_agg+1:end,1:n_agg]
Phi_ss = Phi1[n_agg+1:end,n_agg+1:end]

# partition Se
Se_yy = Se[1:n_agg,1:n_agg]
Se_ys = Se[1:n_agg,n_agg+1:end]
Se_sy = Se[n_agg+1:end,1:n_agg]
Se_ss = Se[n_agg+1:end,n_agg+1:end]

# Initialization
s[1,:]   = s_obs[1,:];
P[1,:,:] = H_t[:,:,1];

# Kalman Filter Recursion
liki     = zeros(T,1);

# Forward Iterations
for t=2:T

    # forecast state and yobs
    yprime = Phi_yy*y_obs[t-1,:]+Phi_ys*s[t-1,:];
    sprime = Phi_sy*y_obs[t-1,:]+Phi_ss*s[t-1,:];

    Pprime_yy = Phi_ys*P[t-1,:,:]*Phi_ys' + Se_yy;
    Pprime_sy = Phi_ss*P[t-1,:,:]*Phi_ys' + Se_sy;
    Pprime_ss = Phi_ss*P[t-1,:,:]*Phi_ss' + Se_ss;

    # condition distr of state on y_obs
    sprime_y_obs    = sprime + Pprime_sy*inv(Pprime_yy)*(y_obs[t,:]-yprime);
    Pprime_ss_y_obs = Pprime_ss - Pprime_sy*inv(Pprime_yy)*Pprime_sy';

    # forecast observation
    yobs_pred = yprime
    sobs_pred = sprime_y_obs;
    v_y       = y_obs[t,:] - yobs_pred;
    v_s       = s_obs[t,:] - sobs_pred;
    F_yy      = Pprime_yy;
    F_ss_yy   = Pprime_ss_y_obs + H_t[:,:,t]

    # if det(F) <= 0.0
    #     U,Sig,V = svd(F)
    #     detF = prod(Sig[Sig.>1e-12])
    # else
    #     detF = det(F)
    # end
    liki[t]  = ( -0.5*n_agg*log(2*pi) - 0.5*log(det(F_yy)) - 0.5*v_y'*inv(F_yy)*v_y
                 -0.5*n_s*log(2*pi) - 0.5*log(det(F_ss_yy)) - 0.5*v_s'*inv(F_ss_yy)*v_s)

    # updating step
    kgain      = Pprime_ss_y_obs*inv(F_ss_yy);
    s[t,:]     = sprime_y_obs + kgain*v_s;
    P[t,:,:]   = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;

end

return s, liki

end


#-----------------------------------------------------------

#-----------------------------------------------------------

function KFwithX(H_t,Se,Phi1,y,n_agg,densdraws,employdraws,knots,xgrid,coef_lambda,coef_mean)

(T,n) = size(y)
n_s   = n-n_agg
s     = zeros(T,n_s);
P     = zeros(T,n_s,n_s);

y_obs = y[:,1:n_agg]
s_obs = y[:,n_agg+1:end]

# partition PHI1
Phi_yy = Phi1[1:n_agg,1:n_agg]
Phi_ys = Phi1[1:n_agg,n_agg+1:end]
Phi_sy = Phi1[n_agg+1:end,1:n_agg]
Phi_ss = Phi1[n_agg+1:end,n_agg+1:end]

# partition Se
Se_yy = Se[1:n_agg,1:n_agg]
Se_ys = Se[1:n_agg,n_agg+1:end]
Se_sy = Se[n_agg+1:end,1:n_agg]
Se_ss = Se[n_agg+1:end,n_agg+1:end]

# Initialization
s[1,:]   = s_obs[1,:];
P[1,:,:] = H_t[:,:,1];

# Kalman Filter Recursion
lh_ss     = zeros(T,1);
lh_X      = zeros(T,1);

# Forward Iterations
for t=2:T

    # forecast state and yobs
    yprime = Phi_yy*y_obs[t-1,:]+Phi_ys*s[t-1,:];
    sprime = Phi_sy*y_obs[t-1,:]+Phi_ss*s[t-1,:];

    Pprime_yy = Phi_ys*P[t-1,:,:]*Phi_ys' + Se_yy;
    Pprime_sy = Phi_ss*P[t-1,:,:]*Phi_ys' + Se_sy;
    Pprime_ss = Phi_ss*P[t-1,:,:]*Phi_ss' + Se_ss;

    # condition distr of state on y_obs
    sprime_y_obs    = sprime + Pprime_sy*inv(Pprime_yy)*(y_obs[t,:]-yprime);
    Pprime_ss_y_obs = Pprime_ss - Pprime_sy*inv(Pprime_yy)*Pprime_sy';

    # forecast observation
    yobs_pred = yprime
    sobs_pred = sprime_y_obs;
    v_y       = y_obs[t,:] - yobs_pred;
    v_s       = s_obs[t,:] - sobs_pred;
    F_yy      = Pprime_yy;
    F_ss_yy   = Pprime_ss_y_obs + H_t[:,:,t]

    # if det(F) <= 0.0
    #     U,Sig,V = svd(F)
    #     detF = prod(Sig[Sig.>1e-12])
    # else
    #     detF = det(F)
    # end
    lh_ss[t]  = ( -0.5*n_agg*log(2*pi) - 0.5*log(det(F_yy)) - 0.5*v_y'*inv(F_yy)*v_y
                 -0.5*n_s*log(2*pi) - 0.5*log(det(F_ss_yy)) - 0.5*v_s'*inv(F_ss_yy)*v_s)

    # reconstruct likelihood function for cross-sectional observations
    densdraws_t     = densdraws[1:N,t]
    employdraws_t   = employdraws[1:N,t]
    selecteddraws_t = densdraws_t[employdraws_t.==1]
    coef_t          = coefRecover(s_obs[t,:]', coef_lambda, coef_mean)
    lh_X[t]         = -length(selecteddraws_t)*logspline_obj(selecteddraws_t, coef_t', knots, minimum(xgrid), maximum(xgrid))

    # updating step
    kgain      = Pprime_ss_y_obs*inv(F_ss_yy);
    s[t,:]     = sprime_y_obs + kgain*v_s;
    P[t,:,:]   = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;

end

return s, lh_ss, lh_X

end
